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SUMMARY 


This  paper  concerns  principles  for  the  selection  of  variational  discretiza¬ 
tion  methods.  The  effectiveness  of  a  method  is  analyzed  in  terms  of  the  key 
concepts  of  approximability  and  stability  which  are  related  to  the  goals  of  the 
computation  through  a  selected  norm  on  the  error.  The  importance  of  robustness 
is  emphasized  and  its  influence  on  the  selection  of  methods  is  analyzed.  These 
concepts  are  elaborated  on  through  a  series  of  illustrative  examples  and  results 
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.1  INTRODUCTION 


The  goal  of  engineering  computations  is  to  obtain  quantitative  information 
about  engineering  problems.  This  goal  is  usually  achieved  by  the  app roximate  solution 
of  a  mathematically  formulated  problem.  Although  a  relevant  mathematical  for¬ 
mulation  of  the  problem  and  its  approximate  solution  are  closely  related  (see 
e.g.  [1]*[2]),  here  we  shall  suppose  that  a  mathematical  formulation  has  already 
been  determined  and  is  amenable  to  an  approximate  treatment.  We  shall  discuss 
a  broad  class  of  approaches  based  on  variational  methods  of  discretization  which 
allow  one  to  find  the  approximate  solution  within  a  desired  range  of  accuracy. 

Let  H  denote  the  linear  space  of  possible  solutions  and  u€H  the  exact 
solution  of  the  problem.  A  (linear,  consistent)  variational  method  of  discretiza¬ 
tion  consists  of  a  finite  dimensional  linear  subspace  SCH  called  the  trial 
space  in  which  the  approximate  solution  is  sought,  a  test  space  V  (of  the  same 
dimension  as  the  trial  space  S  ) ,  and  a  bilinear  form  B(u,v)  defined  on  H  *  V  . 
The  approximate  solution,  denoted  by  Pu  ,  is  then  determined  by  the  conditions 

(  -1.1a)  Pu  €S  , 

(  .1.1b)  B(Pu,v)  =  B(u,v)  ,  for  all  v£V  . 

In  order  that  Pu  be  computable,  the  following  two  conditions  must  be 
satisfied : 

(  .1.2a)  for  any  v€V  ,  B(u,v)  is  computable  from  the  data  of 

the  problem  (without  knowing  u  )  ; 

(  .1.2b)  for  any  s€S  there  is  some  v€V  such  that  B(s,v)  i*  0  . 
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It  follows  that  (  .1.1)  leads  to  a  system  of  linear  equations  which  is  uniquely 
solvable.  It  is  obvious  that  Pu  =  u  for  any  u€S  .  The  approx¬ 
imate  solution  Pu  obviously  depends  on  the  selection  of  S  ,  V  ,  and  B  . 

The  acceptability  of  the  approximate  solutionis  stated  in  terms  of  a  norm 
| I'll  of  the  difference  of  u  and  Pu  ,  i.e.,  we  accept  Pu  if 

(.1.3)  || Pu-u | |  £  t | | u | |  , 

where  t  is  an  a-priori  given  tolerance.  (An  absolute  error  criterion  or  other 
variant  is  equally  possible.)  Thus,  given  ||‘||  and  x  ,  the  goal  is  to  select 
S,  V,  and  B  so  that  (  .1.3)  is  achieved  in  a  most  effective  way.  (We  do  not 
give  here  an  exact  meaning  tj  the  word  "effective".) 

For  each  mathematical  problem  there  exists  a  wide  variety  of  possible 
vanattcnal  methods  of  discretisation.  In  this  paper  we  shall  discuss  properties 
of  these  methods  which  enable  us  to  distinguish  among  them  and  which  aid  in  the 
selection  or  design  of  a  method  which  is  effective  in  achieving  the  given  goals 
of  the  computation.  In  Section  .2  some  general  considerations  are  discussed.  The 
remainder  of  the  paper  is  devoted  to  specific  illustrative  results.  We  conclude 
this  introduction  with  a  very  simple  example  in  terms  of  which  some  of  the  main  ideas 
will  be  explained. 

Let  us  consider  a  longitudinally  loaded  bar  on  an  elastic  support.  For 
0  <  x  <  l  denote  by  u(x)  and  a(x)  the  longitudinal  displacement  and  normal 
stress,  respectively.  A  classical  formulation  consists  of  the  boundary  value 
problem 


(  .1.4a) 


E(x)u*  (x)  ■*  o(x) 
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(  .1.4b)  -(F(x)a(x))'  +  c(x)u(x)  =  p(x)  , 

(  .1.4c)  u(0)  =  0  ,  u ( SI)  =  0  . 


Here  E(x)  denotes  the  modulus  of  elasticity,  F(x)  the  crossection,  c(x)  the 
spring  constant  of  the  elastic  support,  and  p(x)  the  longitudinal  load.  We 
can  cast  this  problem  in  a  variational  form  in  various  ways.  For  example,  define 
the  bilinear  form 

% 

(  .1.5)  B..(u,a;v,p)  =  (Eu'p-op+Fav'+cuv)dx  . 

1  J  0 


Then 


and  so  B^(u,a;v,p)  is  computable  without  explicitly  knowing  the  exact  solution 

°1  0 

(u,  a).  Note  that  B^(u,a;v,p)  is  defined  for  all  (u,a)€H  x  H  =  H  and  all 
(v,p)  €H  . 

There  are  many  other  bilinear  forms  which  could  be  used  in  a  variational 
formulation  of  (  .1.4).  For  example  let 


(  .1.7) 


B»(u,a;v,p)  =  (up - +  Fav  +  cuv)dx 

2  0  E 


(u,a;v,p)  =  | 


pv  dx 
0 


for  (u,a)  ,  (v,p) €  H  (with  H  defined  as  above).  Integrating  by  parts  in 
(  .1.7)  we  get 


X 

(  .1.8)  B~(u,o;v,p)  =  [-up’  -  -^r-  -  (Fo)'v  +  cuv)dx 

J  Jn  E 


B3(u, 


x 

o;v,p)  = 

;0 


pv  dx 


Here  we  assume  u  ,  v€  H  and 


x 

f.pCH1  =  {v|j  [v2  +  (v')2]dx  <  <*>} 


In  both  cases  the  bilinear  form  is  obviously  computable  from  data.  For  a  final 
example  set 


(  .1.9) 


x 

i,v)  =  j  (EFu’v' 
1  0 


+  cuv)dx 


By  eliminating  o  from  (  .1.4)  we  see  that 


(u»v)  *  | pv  dx 


for  u,v€n  .  This  is  the  usual  form  used  in  displacement  finite  element  methods 
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Of  course  many  other  variational  formulations  of  (  .1.4)  are  possible  (in  fact 
an  infinite  number). 

To  complete  the  specification  of  a  discretization  method  we  must  select  in 
addition  to  the  bilinear  form,  the  trial  and  test  spaces  S  and  V  .  For 
example,  in  the  case  of  we  may  select  any  finite  dimensional  subspaces 

S  and  V  of  H  which  are  of  the  seme  dimension  and  satisfy  (  .1.2b). 


Let  us  now  define  some  of  the  norms  which  we  will  consider.  For  k  a 

,  [k]  d\i 

nonnegatxve  integer  set  u  =  — r  and  let 

dx 


Y  ess  sup  |u»l  (x)  | 
j  =0  0<x<£ 


We  also  define  analagous  norms  for  k  a  negative  integer.  For  such  k  define 


[k]  . 

v  =  u  by 


vf2m]dx=0. 


u (x)  =  v^  ^(x)  , 
0  <  m  £  , 


choosing  the  constants  of  integration  so  that 
v(2nri-l](0)  =  v[2mfl]w  _  0  >  0  £  m  <  =~-  . 


Then  we  will  write 
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These  negatively  indexed  norms  emphasize  the  effect  of  oscillations  of  u  less 
than  do  the  positively  indexed  norms.  Analogous  norms  can  be  defined  when  u 
depends  on  two  variables. 


.2  GENERAL  CONSIDERATIONS 

To  achieve  the  acceptance  criterion  (  .1.3)  it  is  certainly  necessary  that 

(  .2.1)  Z(u,S)  =  inf ] | u-s | |  <  T 

s€S 

The  quantity  Z(u,S)  measures  the  error  in  best  possible  approximation  of  u  by 
elements  of  S  with  respect  to  the  chosen  norm  ||-||  ,  i.e.,  the  approximability 
of  u  by  S  . 

The  choice  of  *'  is  clearly  essential  to  the  effectiveness  of  the  dis¬ 
cretization  method.  The  solution  u  is  unknown  a-priori,  and  often  only  the 
information  that  u€H  is  available.  In  such  cases  S  has  to  be  selected  so 
that  every  element  in  H  can  be  approximated  well.  More  Information  about  u 
allows  more  effective  choice  of  S  .  Such  information  can  be  achieved  through 
a  learning  process  during  the  computation  and  thus  S  can  be  selected  adaptively. 
See,  e.g.,  [3] ,  [4] . 

That  the  trial  functions  approximate  the  solution  well,  i.e.,  that  the 
magnitude  of  Z(u,S)  is  small,  does  not  alone  insure  that  the  approximate 
solution  pu  is  close  to  the  exact  solution  u  .  Therefore  it  is  reasonable  to  ask 
that  the  method  be  quasioptimal .  This  means  that 

(  .2.2)  |  |  u— Pu  |  |  £  KZ(u,S)  =  K  inf  ]  |  u-s  |  |  ,  for  all  u€H  , 

s€5 


where  K  is  a  constant  which  is  not  too  large.  The  smallest  value  of  K  for 


which  (  .2.2)  holds  is  called  the  quasioptimality  constant. 


Condition  (  .2.2)  is  equivalent  to  another  condition  called  the  stability 
condition.  This  states  that 


(  .2.3)  |  j Pu |  |  £  K  |  |u|  j  ,  for  all  u€H 

* 

The  smallest  value  of  K  for  which  (  .2.3)  holds  is  called 

the  stability  constant.  To  see  that  quasioptimality  and  stability  are  equivalent, 
assume  that  (  .2.3)  holds.  Now,  if  s€S  ,  then 

I  I u— Pu | |  =  | | (u-s)-P(u-s) | )  £  j |u-s| ]  +  | |P(u-s) | | 

<  I |u-s I |+K  I |u-s| I  £  (K  +1) I |u-s I  I 

(Here  we  used  the  fact  that  Ps  =  s  ,  as  mentioned  earlier.)  Thus  (  .2.2)  holds 
* 

with  K  <  K  +  1  .  On  the  other  hand  assuming  (  .2.2)  we  have  that 

IIPull  £  I  I Pu-u I  I  +  ||u||  £  KZ(u, S)  +  ||u||  £  (K+l)||u||  , 

* 

and  so  (  .2.3)  holds  with  K  <  K  +  1  .  The  importance  of  (  .2.3)  is  that  it  is 
often  easier  to  verify  then  (  .2.2). 

Note  that  while  approximability  is  affected  only  by  the  choice  of  the  trial 
space  S  ,  stability  (or  quasioptimality)  depends  on  the  interplay  between  B  , 

H  ,  S  ,  and  V  .  Because  the  test  functions  are  not  needed  for  approximation 
purposes,  the  main  goal  in  the  selection  of  V  is  to  achieve  stability  with  the 
smallest  possible  constant  K  .  Let  us  remark  that  for  certain  bilinear  forms 
and  certain  norms,  the  choice  V  *  s  leads  to  stability  constant  1  .  In 


j 
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such  cases,  the  performance  of  the  method  depends  solely  on  the  selection  of  S  . 

Note  further  that  both  approximability  and  stability  depend  heavily  on  the 
norm  under  consideration.  Changing  the  norm  can  violate  quasioptimality  although 
the  computational  algorithm  remains  the  same.  Because  of  this,  the  method  must 
be  investigated  in  close  relation  to  the  given  acceptance  criterion. 

Although  approximability  and  stability  are  essential  and  of  primary  interest 
for  the  method,  there  are  other  important  features  to  be  considered  in  the  rational 
selection  of  discretization  procedures. 

Robustness .  The  bilinear  form  B  and  the  solution  u  may  depend  on 
various  parameters,  e.g.,  in  the  above  example  of  the  bar  problem,  E  ,  F  ,  and 
c  may  play  a  significant  role.  Both  approximability  and  stability  will  depend 
on  such  parameters.  A  method  is  called  robust  when  its  performance  is  relatively 
uninfluenced  by  the  variation  of  the  parameters  within  a  large  range. 

A-posteriori  estimates  and  adaptive  approaches.  A  typical  acceptance 
criterion,  as  mentioned  above, is 

I |u-Pu| I  <  t| |u| I  , 
where  r  is  a  given  tolerance.  Although  we  have 

j |u-Pu| J  £  KZ(u,S)  , 

this  estimate  may  have  no  direct  practical  importance.  In  the  first  place  we 
will  in  general  not  know  precise  values  for  the  quasioptimality  constant  K  or 
for  Z(u,S)  .  Moreover  even  when  these  are  known  the  resulting  estimate  may  be 
very  pessimistic.  The  reason  is  that  the  quasioptimality  constant  K  is  based 
on  the  worst  case,  (since  (  .2.2)  must  hold  for  all  uCH),  while  the  true  solu¬ 
tion  may  have  special  properties  unknown  to  us.  The  only  general  ways  to  implement 
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the  acceptance  criterion  reliably  are  based  on  a-posteriori  analysis  of  the  ap¬ 
proximate  solution  l’u  .  Thus  a  computable  error  estimator  e  is  introduced, 
which  depends  solely  on  input  data  and  Pu  and  satisfies 

e  -v  I  j  u-Pu  |  | 


in  the  sense  that 


and 


C^c  £  | [ u-Pu | |  £  C2 e 

A  A 


0 


e _ 

|  u-Pu ] j 


->  1 


as 


c  -*•  0 


This  can  be  achieved  (see  e.g.  [3],  [4],  [5],  [6]),  but  not  every  selection  of  H  ,  S  , 
V  ,  B  ,  and  | | ’ | |  allows  for  estimators  with  the  same  effectiveness  and  reliability. 
Feasibility  of  adaptive  selection  of  test  and  trial  spaces  may  also  be  an  important 
feature  to  be  considered  in  the  selection  of  the  form  B  . 


.3  ILLUSTRATIVE  RESULTS 

In  this  section  we  discuss  some  concrete  mathematical  results  illustrating 
the  ideas  introduced  above. 

.3.1  Approximal ility 

First  we  consider  some  questions  related  to  approximation.  In  engineering 
computations  the  solution  we  are  interested  in  approximating  usually  has  special 
properties.  For  example  it  may  be  smooth  except  for  some  singular  behavior 
in  the  neighborhood  of  a  knc*m  point  such  as  a  crack  tip,  comer,  or  concentrated 
load.  Moreover  the  qualitative  nature  of  such  singularities  is  known. 
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•3.1.1  The  one  dimensional  case 

The  one  dimensional  analog  of  "comer"  behavior  in  2  and  3  dimensions  is 

y 

given  by  functions  u^(x)  =  x  ,  y  a  real  number.  This  function  has  the  property 
that 


(  .3.1) 


k 


l 

1=0 


1 

2£-2k+a 

X 

0 


(d_U) 

V 

dx 


2 

dx  < 


CO 


for  any  integer  k  >  0  and  any  real  number  a  >  2k  -  1  -  2y  .  Suppose  we  are 
interested  in  approximating  in  the  norm  a  function  satisfying  (  .3.1).  It 

can  be  shown  that  there  exists  a  sequence  of  subspaces  of  such  that 

S(n)  has  dimension  n  and  the  satisfy  the  following  approximation 

property:  if  u€H^  is  any  function  satisfying  (  .3.1)  for  a  nonnegative 
integer  k  and  any  real  number  a  such  that  2k  >  a  >_  0  ,  or  if  u€H  ,  then 


(  .3.2) 


Z(u,S(n))  £  C(k  ,a)n  k 


Moreover,  (  .3.2)  exhibits  the  best  rate  of  convergence  achievable  by  any  sub¬ 
spaces  of  dimension  n  .  See  [  7  ].  This  is  a  very  robust  approximation 

property.  In  particular,  all  the  functions  u  ,  with  y  >  -  1/2  can  be  approxi¬ 
mated  with  this  rate.  (For  y  <_  -  1/2,  u  A  so  such  a  result  cannot  apply.) 

In  fact,  since  the  functions  u^  have  additional  properties,  even  better 
approximation  than  indicated  by  (  .3.2),  namely  an  exponential  rate  of  convergence, 
may  be  obtained  for  them  by  another  choice  of  spaces.  Thus  there  exists  a 
sequence  of  subspaces  of  dimension  n  such  that  if  y  >  -  1/2  ,  then 


Z(u^,S(n))  £  Ce_S,/" 


1 


(  .3.3) 
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for  some  B  >  0  .  See  [  8  ]  for  details. 

The  two  results  quoted  above  relate  to  the  existence  of  a  sequence  of  sub¬ 
spaces  of  with  good  approximation  properties.  We  now  consider  the  quality 

of  approximation  achieved  by  some  concrete  choices  of  the  sequences 
suitable  for  computation.  First  let  be  the  space  of  polynomials  of  degree 

less  than  n  .  Then  any  function  satisfying  (  .3.1)  can  be  approximated  with  the 
error 


(  .3.4) 

Applying 


Z(u,P^)  =  inf|ju-si 


sCP 


(n) 


3.4)  to  the  functions  u 


.  .  -min(k,2k-a) 

£  C(a,k)n  . 

,  y  >  -  1/2  we  get  the  estimate 


(  .3.5)  Z(uY,P(n))  <  C(y,e)n'(1+2Y)+e 

with  e  >  0  arbitrarily  small.  It  can  also  be  shown  that  the  estimate  (  .3.5) 
is  essentially  the  best  possible  one. 

Let  us  now  select  ,  the  space  of  all  piecewise  polynomials  of 

degree  less  than  p  on  a  quasi-uniform  partition  of  [0,1]  into  n  elements. 
This  space  has  dimensions  roughly  proportional  to  n  .  The  results  analagous 
to  (  .3.4)  and  .3.5)  in  this  case  are 


(  .3.6) 


Z(u,S^n^)  <  C(k,a,p)n 


-min(p ,k  -  q) 


and 

,  .  -min(p,  -^*-)+e 

(  .3.7)  Z(u  ,SltU)  <  C(Y,e)n 

Y  P  - 


and  these  rates  are  essentially  unimprovable.  Comparing  (  .3.4)  and  (  .3.6)  we 


see  that  for  functions  u  only  assumed  to  satisfy  (  .3.1)  the  rate  of  approxima¬ 
tion  achieved  by  the  polynomials  is  certainly  not  worse  and  may  be  better  than 
that  achieved  by  the  piecewise  polynomials.  For  the  functions  ,  (  .3.5) 

and  (  .3.7)  show  that  the  rate  achieved  by  the  polynomials  is  at  least  twice 
that  achieved  by  the  piecewise  polynomials.  For  more  see  [9]. 

The  estimate  (  ,3.6)  is  in  essence  the  classical  estimate 

Z(u,s^n))  <  C(k,p)n_k|  j  u |  |  k 

when  p  >  k  and  a  =  0  .  See,  e.g.,  [10  ]  •  The  question  arises  whether  under 

these  conditions  an  expression  for  C(k,p)  can  be  given  which  explicitly 

characterizes  the  behavior  with  respect  to  p  .  In  [11]  such  an  expression  is 

given  in  both  the  one  and  two  dimensional  cases  (and  for  approximation  in  H1) . 

A  — k+£  ~ 

It  is  shown  there  that  on  the  right  hand  side  we  can  have  C(k)n  with  C 

independent  of  p  . 

Neither  the  piecewise  polynomial  spaces  nor  the  polynomial  spaces  achieve 

the  optimal  rate  of  convergence  characterized  by  (  .3.2).  For  example  for 

y  >  -  1/2  sufficiently  small,  the  function  u  is  not  approximated  at  the 
-k 

optimal  rate  of  n  for  either  of  these  cases.  Such  a  rate  can  be  achieved 
in  the  first  case  by  a  proper  refinement  of  the  mesh  in  a  neighborhood  of  the 
origin  and  in  the  second  case  by  changing  the  polynomials  to  some  other  system 
of  functions  (see  [11],  [12]).  The  importance  of  this  observation  is  that 
for  engineering  computations  it  appears  likely  that  approximation  spaces  can  be 
created  which  yield  a  rate  of  convergence  which  is  better  than  polynomial,  and  is 
probably  exponentional . 

Thus  far  we  have  considered  approximation  in  the  H°  norm.  Similar  results 
are  available  for  all  the  H  norms,  l  both  positive  and  negative.  An  interest¬ 
ing  fact  is  that  when  i  decreases  the  rate  of  convergence  furnished  by  either 


13 


p<n>  or  S(n)  increases  linearly  with  l  .  For  more,  e.g.,  see  [13]. 

P 

.3.1.2  The  multidimensional  case 

So  far  we  have  discussed  only  the  one  dimensional  case.  Analogous  results 
exist  in  more  than  one  dimension,  but  these  are  far  from  complete.  We  will  not 
go  into  details  here,  but  refer  the  reader,  e.g.,  to  [9  ],  f  H] ,  [14]. 

.3.1.3  The  h-,  p- ,  and  h-p-versions  of  the  finite  element  method 

As  was  said  above,  there  are  important  cases  when  selecting  the  same  trial 
and  test  spaces  leads  to  a  stability  constant  of  1;  and  so  approximability  by 
the  trial  space  determines  the  performance  of  the  method.  The  classical  finite 
element  method  uses  piecewise  polynomials  of  fixed  degree  p  on  meshes  which 
are  refined  to  achieve  accuracy.  Because  the  size  of  the  elements  is  usually 
denoted  by  h  ,  this  method  is  called  the  h-version  and  the  approximability 
properties  (  .3.6)  and  (  .3.7)  for  such  spaces  over  quasi-uniform  meshes  are  used. 
The  p-version  achieves  accuracy  by  fixing  the  mesh  and  increasing  the  degree  p 
of  the  polynomial.  In  this  case  the  approximation  results  (  .3.4)  and  (  .3.5) 
are  applicable.  The  p-version  has  been  implemented  in  the  program  COMET  X  . 

We  refer  the  reader  to  [9],  [15]  and  references  therein  for  detailed  information. 
Finally  the  h-p-version  combines  both  of  these  approaches.  The  exponential 
convergence  rate  given  in  (  .3.3)  can  be  realized  in  the  h-p-version. 

.3.2  Finite  Element  Methods 

We  turn  now  to  a  discussion  of  finite  element  methods.  As  discussed  in 
Chapter  .2,  the  quality  of  approximation  yielded  by  such  a  method  is  assured  by 
stability  in  conjunction  with  approximability.  The  stability  of  a  method  depends 


on  the  interplay  between  the  spaces  S  and  V  ,  the  bilinear  form  B  ,  and  the 
norm  ||*||  •  This  Is  illustrated  in  the  first  example. 
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.3.2.1  An  example  Illustrating  the  role  of  the  trial  and  test  spaces  In  stability 
First  we  consider  a  one  dimensional  problem  with  the  simplest  possible  bi¬ 
linear  form.  Setting  l  =  1  ,  EF  =  1  ,  and  c  =  0  in  (  .1.9)  we  get  the  form 


(  .3.8) 


B(u,v) 


1 

u ' v ' dx  . 

0 


°1 

The  solution  u£H  satisfies  B(u,v) 
two  point  boundary  value  problem  is 


1 

pvdx  for  all 

0 


°1 

v£H 


and  the  related 


-u"  =  p 


u(0)  =  u(l)  =  0 


For  discretization  we  define  spaces  of  smooth  splines.  Let  A  =  {0  =  x^  <  x^  < . . . 
<  =  1}  be  a  mesh  of  [0,1]  and  set  h^  *  x^  -  •  For  y  1  1  ,  the  mesh 

is  called  y-quasiuniform  if 


max  t—  <  y 
,  ,  h .  — 

i.j  J 


A  weaker  restriction  is  that  the  mesh  be  y-locally  quasiuniform,  i.e.,  that 


max 


i  hi±l 


Given  any  mesh  A  we  define  for  r  =  0,1,2,...  the  smooth  splines 

degree  r  subordinate  to  A  to  be  the  piecewise  polynomials  of  degree 
r-1  continuous  derivatives.  The  space  of  all  such  splines  is  denoted 
In  particular  M^(A)  is  the  space  of  piecewise  constant  functions  and 
is  the  space  of  continuous  piecewise  linear  functions.  We  also  denote 


of 

r  with 
Mr(A)  . 
M*(A) 
by  M*(A) 
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the  space  of  piecewise  linear  functions  with  zero  boundary  values  and  by 
M  (A)  the  space  of  natural  cubic  splines,  i.e.,  MJ(A)  =  {v€MJ(A)|v  =  v"  =  0 
when  x  =  0  or  1). 

We  consider  the  use  of  such  spline  spaces  for  S  and  V  in  conjunction 
with  the  bilinear  form  B  defined  in  (  .3.8).  It  is  possible  to  show  that 
if  S  and  V  are  taken  to  be  spaces  of  smooth  splines  of  degree  r^  and  r2 
respectively  (with  appropriate  boundary  conditions)  and  if  r^  and  r 2  are 
either  both  odd  or  both  even,  then  condition  (  ,1.2b)  is  satisfied.  Hence 
Pu£S  is  uniquely  defined  by  (  .1.1). 

The  most  standard  case  occurs  with  S  =  V  =  M^(A)  .  The  stability  properties 

of  this  method  in  several  norms  are  summarized  in  the  following  theorem. 

Theorem  .3.1.  Let  S  *  V  -  M^(A)  for  an  arbitrary  partition  A  .  Then 

(  .3.9a)  | | Pu-u | |  £  K  inf  | | u-s | |  , 

H2  s€S  H2 

(  .3.9b)  |  |  Pu-u  |  |  <_  K  inf  ||u-s||  .  , 

H  s€S  H 

00  OO 

(  .3.9c)  I  I Pu-u ( [  <  K  inf  ||u-s||  n 

HU  s€S  H 

00  ,  OO 

for  all  u€fi^  ,  with  K  independent  of  A  .  However  for  any  C  >  0  and  any 
mesh  A  there  exists  uCH^  such  that 


_ _ 


16 


|  | Pu— u | !  >  C  inf  | |u-s  j 1 

h9  s€S  H? 

4  l. 


□ 


For  more  see  [16]. 


The 

case  where 

S 

=  M*(A)  , 

V  =  M3(A) 

is  less  familiar  and  more  involved 

Theorem 

.3.2.  Let 

s 

=  M1 (A)  , 

V  =  M3(A) 

.  Then : 

a)  For  an  arbitrary  partition  A 


|pu-u| |  0  £  K  inf  | | u-s | |  Q 
H2  s€S  H2 


with  K  independent  of  A  and  u€H 

b)  For  any  Y  ^  1  there  exists  a  constant  K(y)  such  that  for  all 

°1 

y-quasiuniform  partitions  and  all  u€H  , 


|Pu-u||  .  £  K(y)  inf  ||u-slS 

H2  s€S  H2 


and 


! |pu-u! j  ,  <  K(y)  inf  | | u-s | |  , 


s€S 


c)  However  for  any  C  >  0  and  any  partition  A  ,  there  exists  u€H2 


|  pu-u  |  j  ~  C  inf  j  |  u-s  |  | 


s€S 


2 


I 


such  that 
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d)  If  1  <  y  <  y  =  l+  /3  +/ 3+2/3  =  5.2745...,  there  exists  a  constant 


K(y)  such  that 


and 


I  Pu-u | |  .  <  K(y)  inf  |  | u-s | | 

H2  s€S  H2 


|  Pu-u | |  <  K(y)  inf  | | u-s | |  _ i 


s€S 


for  any  Y-locally  quasiuniform  partition  A  and  u  .  However 
lim  K(y)  =  00  •  I  I 

Y»Y0 

Thus  we  see  that  there  is  a  very  fine  interplay  between  the  trial  and  test 
spaces  even  for  the  simplest  bilinear  form. 

So  far  we  have  analyzed  the  form  (  .3.8).  We  now  consider  the  form 


(  .3.10) 


j. 

B(u,v)  =  |  Eu'v'dx 


where  0  «.  eQ  £  E(x)  <  <  ”  .  The  question  arises  whether  Theorems  .3.1  and 

.3.2  remain  true  as  stated.  It  is  possible  to  show  that  if  E(x)  is  suf¬ 
ficiently  smooth,  then  Theorems  .3.1  and  .3.2  hold  without  change.  The 
requirement  of  smoothness  means  that  K  may  depend  as  well  as  on  e^  and  e^, 
also  on  the  maximum  of  the  first  few  derivatives  of  E  .  It  can  be  shown  that 
(  .3.9a)  holds  with  K  depending  only  on  eg  and  e^  ,  but  (  .3.9c)  is  not  true 
when  no  differentiability  restrictions  are  made  on  E  .  Thus  we  may  say  that 
the  performance  of  that  method  is  more  robust  with  respect  to  the  coefficient  E 

j ) • | |  than  it  ia  in  the  norm  |  |  •  |  |  n 
Hn  K 

y  00 


in  the  norm 
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.3.2.2  The  bilinear  form  and  robustness 

We  continue  to  consider  the  simple  problem  (  .1.4)  but  now  consider  the 
effect  of  the  choice  of  the  bilinear  form  on  the  robustness  of  the  method.  For 
simplicity  assume  that  F  =  1  and  c  =  0  .  Let  E(x)  be  given  satisfying 
0  <  e^  <  E(x)  <  e^  and  consider  the  bilinear  forms  (  .1.5)  and  (  .1.7).  The 
bilinear  form  (  .1.5)  clearly  stems  from  the  system  of  equations 

(  .3.11)  Eu'  =  a  , 

o'  =  p  , 

while  (  .1.7)  comes  from 

(  .3.12)  u* 

a' 

Assume  now  that  we  take 

■  1  0 

(  .3.13)  S  =  V  =  Mx (A)  x  m  (A) 

in  both  cases.  It  is  easy  to  see  that  o  can  be  eliminated  from  the  system  of  linear 
equations  arising  from  (  .1.5)  with  the  choice  of  spaces  given  in  (  .3.13),  and 

.  i 

we  then  get  the  same  method  as  when  (  .3.10)  is  used  with  S  =  V  =  M  (A)  .  The  prop¬ 
erties  of  this  method  were  summarized  in  Theorem  .3.1.  The  form  .(  .1.7)  however 
gives  different  results,  which  we  now  consider  in  detail. 

Letting  (u,a)  *  P(u,a)  we  get 

Theorem  .3.3.  Let  A  be  an  arbitrary  mesh.  Define  P  by  the  bilinear  form 


(  .1.7)  with  S  and  V  defined  by  (  .3.13).  Then  there  exists  a  constant  K 


depending  only  on  and  such  that 


a) 


b) 


I  |u-u|  |  +  ||o-o||  <K  inf  ["||u-xi|  ,  +  ||o-*||  Q]  , 

H2  H2  (X,iJJ)63L  H2  H2J 


||u-u|!  +  |  | o-o  j  |  <  Kp  inf  |  | u-x  1  |  0+  ||o-i|»||  Q  1  , 

fT  H  kX,^€S  H  H 

OO  oo  w  00  oo  J 


c) 


inf 


i|^°(A) 


I  Ml  I  o 

H„ 


d) 


o-a 


<  K 


J3i 


a— 


ij/€Mu(A) 


□ 


The  statements  analogous  to  c)  and  d)  for  the  error  | |u-u| |  are  not  true. 

In  order  to  elaborate  this  point  let  us  introduce  a  further  notation.  For 
°1 

X  €M  (A)  let  x  be  defined  by 

i  n. 

X  (Xj )  =  X(xj)  ,  j  =  0,1,. ..,n, 

(Ex')'  =0  on  »  J  =  l,2,...,n 

Then  we  have 

Theorem  .3.4.  Let  u  be  defined  as  in  Theorem  .3.3  (i.e.,  using  the  form 
(  .1.7)  and  spaces  of  (  .3.13)).  Then 
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|u-u||  o  £  K(eQ,e  ,V(E))  inf  [ | | u-x I |  Q  +  I |u~x| I  Q] 

Hoo  y^M^A)  H”  H 


and 


max | u(x . )-u(x  ) j  £  K(e0,el  ,V(E))  inf  J  j u— x I  I  q  , 
j  xe^iA)  H- 

with  K  depending  on  e^,  ej  ♦  an<^  the  variation  V(E)  of  the  function  E  .  Q 

We  remark  that  this  theorem  is  not  valid  when  the  dependence  of  K  on 

V(E)  is  supressed.  Moreover,  while  the  term  inf||u~x||  n  in  the  second 

H 

oo 

estimate  is  necessary,  it  is  usually  smaller  than  the  first  term.  Comparing 
Theorem  .3.4  with  the  previous  results  we  see  that  the  form  (  .1.7)  is  much 
more  robust  than  (  .1.5)  with  respect  to  all  the  norms  we  have  considered 

except  the  norm  and  should  be  preferred  in  most  situations. 

Let  us  comment  on  the  system  of  linear  equations  which  the  approximate 

solution  lead  to  when  o  is  eliminated  in  either  of  the  two  methods  discussed 

•  i 

in  this  section.  In  both  cases  Pu€M  (A)  is  defined  by  a  system  of  the  form 


1 


EA(Pu),v’dx  = 


pvdx,  for  all  v€M^(A) 


0  0 
When  the  form  (  .1.5)  is  used 


EA  h, 


while  when  (  .1.7)  is  used 


xi 


E(x)dx  on  (x1_1,xi)  , 


i-1 


xi  -i-1 

\  E(10dxJ  °n  (xi-l’Xi) 


'i-1 


i”  l,2,...,n  .  Thus  in  the  former  case  E  is  replaced  by  its  piecewise 


average  and  In  the  latter  by  its  piecewise  harmonic  average.  The  above  stated 
results  show  that  the  usual  finite  element  method  does  not  have  as  good  stability 
properties  when  E(x)  changes  significantly  over  an  element,  and  so  should  not 
be  used.  The  change  can  be  measured  by  the  ratio  of  the  average  to  the  harmonic 
average . 

.3.2.3  Changing  the  dependent  variable  to  improve  approximability 

Now  we  turn  to  the  analysis  of  the  form  defined  in  (  .1.8)  where  for 
simplicity  we  take  F  =  1  and  c  =  0  .  If  we  choose  S  =  V  =  M°(A)  x  M1(A)  and  set 
(u,a)  =  P(u,a)  ,  it  is  easy  to  prove  that 


(  .3.14)  | |u-u| |  +  | | ct-ct I |  <  C 

H  H 

with  C  independent  of  A  but  depending  on  E  .  Now  for  any  gCH1 
consider  the  variational  formulation 

l  i 

(  .3.15)  B„(u  ,a  ;v,o)  =  [  (p-g')vdx  -  f  ^  dx 

8  8  )0  L 


instead  of  the  method  just  considered: 


Bj(u,o  ;v,p) 


l 

pvdx 

0 


The  new  variables  are  related  to  the  old  by 
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Thus  we  may  compute  a  and  then  take  a  =  a  +  g  for  the  stress  component  of 

8  8 

the  approximate  solution.  Because  the  same  bilinear  form  occurs  in  both  these 
approaches,  the  stability  is  unchanged  and  we  get 

Theorem  .3.5  For  the  method  associated  with  (  .3.15) 


where  the  C  is  the  same  constant  which  appears  in  (  .3.14)  (and  therefore  is 

independent  of  g  ) .  Moreover  a  -  a  =  o  -a  .  (Hj 

8  8 

Now  we  note  that  the  propoer  choice  of  g  can  increase  the  smoothness  of 
0  ,  thereby  increasing  its  approximability  and  so  improving  the  accuracy  of  the 

O 

method.  In  this  simple  one  dimensional  case  the  best  choice  is  g  =  J  pdx  so 
0g  is  constant.  (In  some  situations  it  is  as  easy  or  easier  for  the  user  to 
input  g  as  p  .)  In  this  case  the  last  term  in  (  .3.16)  will  disappear 
and  o  will  exactly  equal  0  .  A  similar  idea  may  be  fruitfully  applied 
to  related  mixed  methods  in  more  than  one  dimension. 

.3.2.4  A  robust  method  for  a  parameter-dependent  problem 

In  Subsection  .3.2.2  we  discussed  a  simple  case  in  which  changing  the  bi¬ 
linear  form  increased  significantly  the  robustness  of  the  method.  We  now  discuss 
another  example  in  which  a  parameter  enters  in  a  direct  fashion.  The  problem 
to  be  considered  models  the  deflection  of  a  beam  allowing  for  the  effect  of  shear 
stress.  In  the  simplest  case  this  model  can  be  described  by  the  system  of 
equations 


(  .3.17) 


+  d~2Ud-u£)  =  0  ,  0  <  x  <  1  , 

d  2(Vo,-)*  =  g  ,  0  <  x  <  1  , 
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with  the  boundary  conditions 


i>d  (0)  =  4>d(D  =  ..d(0)  -  ud(i)  -  o  . 


Physically  <J>d(x)  represents  the  displacement,  wd(x)  the  rotation  of  the 
cross  section,  and  g(x)  the  transverse  load.  The  solution  depends  on  the 
beam  thickness  d  .  We  associate  to  the  problem  (  .3.17)  the  bilinear  form 


x 

Bd(4>  ,v)  =  [  <f> *  tji'+d  2  ( 4>—  o ’  )  (4j-v')]dx 


Let  S  =  V  =  M1(A)  x  M^(A)  .  Then  we  have 


(  .3.18) 


Bd(<t,d,tVf',V*  =  j  gvdx 
0 


and  so  the  bilinear  form  is  computable.  Denoting  P($>  a^)  by  (<J>d,wd)  we  get 


Theorem  .3.6. 


>d"*d*l  1+  ilv“d>>u  1  -c(d)r  i"! Mlv^li  +  l'Vp|lHi 

H  H  (x»TV  €S  L  H  H  __ 


The  constant  C(d)  is  independent  of  4>,  u>€H  ,  but  C(d)  -*•  as  d  -*•  0  .  | _ j 

A  corollary  of  this  theorem  is 

Theorem  .3.7.  Let  any  nonzero  load  g  be  given  and  let  0  <  a  <  1  be 
arbitrary.  Then  for  any  partition  A  there  exists  a  value  of  d  depending  on 
A  such  that 


24 


I 1 1  i  -  °l 1*1 1  i 

H  H 


lwd-(i'd '  1  1  1  °!  I^d !  1  1 

H  H 


□ 


Theorem  .3.7  shows  that  for  small  d  the  method  based  on  (  .3.18)  is  virtually 
useless. 

We  now  associate  to  our  problem  another  bilinear  form,  in  which  we  introduce 
a  new  variable  £  ,  representing  the  shear  stress. 


(  .3.19)  bdU,us,S;'|J,v,n)  = 


[(fc'^'+^^-v’)  +  n(<t>-a>')  -  d  £n]dx 


0 


-2 

The  functions  <f>  ,  w  and  £ ,  =  d  (<{>.-w')  satisfy 


(  .3.20) 


Bd^<*,d’Wd’?d’^’v,r^  = 


gvdx 


for  any  i^H1  ,  vCH1  ,  n  €H°  .  Select  now  S  =  V  =  MX(A)  x  MX(A)  *  M°(A)  , 
and  let  (4>5(*)»0  *  P(<t>d>wd ’^d^  *  The  robustness  of  the  new  method  with  respect 

to  the  parameter  d  is  evidenced  by  the  following  result. 


Theorem  .3.8.  For  the  method  associated  with  (  .3.20) 


IVVl  1  +  ll“d-“dllui  +  llv^d'i  0  i 

H  H  H 


lnLc[,,Vx|,Hi  +  NvpNi  +  Ncd-MI  J 
(x,p,«s  H  H  H 


with  C  independent  of  A  and  d 


□ 
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When  g  in  (  .3.17)  is  smooth,  then  (Ji^,  and  are  smooth  also, 

and  may  be  well  approximated  independently  of  d  .  It  follows  that  computations 
based  on  (  .3.19)  and  (  .3.20)  give  very  good  results  while  we  have  seen  that 
computations  based  on  (  .3.18)  yield  extremely  poor  results  for  small  d  .  This 
difference  in  the  robustness  of  the  two  methods  with  respect  to  d  is  very 
striking  in  practice.  It  is  also  worth  noting  that  the  additional  variable 

may  be  eliminated  from  (  .3.20).  The  resulting  method  is  identical  with  the 
method  based  on  (  .3.18)  except  that  the  integrals  are  calculated  by  the  com¬ 
posite  midpoint  rule.  By  employing  this  reduced  integration  implementation,  the 
mixed  method  entails  no  extra  expense  whatever.  For  more,  see  [17]. 

.3.2.5  Robust  methods  in  two  dimensions 

So  far  we  have  discussed  various  ideas  concerning  the  proper  selection  of 
S,  V,  and  B  in  the  context  of  simple  one  dimensional  examples.  Analogous  ideas 
can  be  used  in  more  dimensions  also,  although  much  less  is  known  at  present. 
Nevertheless  we  will  briefly  consider  some  examples.  Consider  the  problem 

3  3u  ,  3  3u  c 

3x  3  3x  3y  3  3y 

on 

ft  =  { (x,y) |  | x |  <  1  ,  | y |  <  1 } 

with  the  condition  u  *  0  on  3ft  .  Assume  that  a  *  a^  for  x  <  0 
and  a  =  for  x  >  0  ,  where  Sq  and  a^  are  distinct  positive  constants. 

Having  selected  a  triangulation  T  (with  minimal  angle  condition),  the  usual 
method  employs  the  bilinear  form 


x_. 
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(.3.21,  B<u,v>  -  fa(.£  ♦  .  *  ^dxdy 

and  equal  trial  and  test  spaces  consisting  of  functions  which  are  continuous, 
linear  on  every  triangle,  and  zero  on  3C2  .  If  the  interface  x  =  0  does  not 
coincide  with  the  boundaries  of  the  triangles,  then  the  solution,  which  is  not 
smooth  at  the  interface,  will  be  approximated  less  accurately  than  for  a  problem 
with  a  smooth  coefficient.  We  shall  show  how  to  proceed  in  a  slightly  different 
fashion  which  will  avoid  this  problem. 

We  select  for  V  the  space  of  continuous  piecewise  linear  functions.  Thus 
the  restriction  of  a  test  function  to  a  triangle  is  a  linear  combination  of 
the  three  functions  1  ,  x  ,  and  y  and  each  test  function  is  continuous  at  the 
nodes.  The  trial  functions  are  taken  to  restrict  on  each  triangle  to  a  linear 

f  X  "I 

combination  of  the  functions  1  ,  J  —  dt  ,  and  y  and  to  be  continuous  at  the 

0  a 

nodes.  The  trial  functions,  unlike  the  test  functions,  need  not  be  continuous 
on  element  boundaries  except  at  the  nodes  (and  so  are  "nonconforming") .  Now 
interpret  (  .3.21)  as  a  sum  of  integrals  over  the  individual  triangles  and  re¬ 
place  the  norm  |  |  '  |  I  i  by  the  norm  |  |  ’  |  |  defined  as  the  square  root  of 

Hj  H‘(T) 

a  sum  of  integrals  over  the  triangles.  This  is  a  common  approach  in  nonconforming 
finite  element  methods. 

Theorem  ,3.9j  For  this  method 

| u-Pu | |  .  £  C  inf  | ( u-X  j |  . 

HV)  X€S  hV) 

with  C  independent  of  T  and  u  .  The  value  of  the  unusual  trial  space  S 
used  here  is  that  while  stability  still  holds  (as  stated  in  Theorem  .3.9),  these 
trial  functions  mimic  the  behavior  of  the  solution  u  and  thus  greatly  improve 


the  approximability .  That  is,  Z(u,S)  is  generally  much  smaller  for  this  choice 
of  S  than  if  S  is  taken  equal  to  V  (the  usual  choice).  The  resulting  method 


therefore  gives  superior  results.  An  important  observation  to  be  made  here  is 
that  the  difficulties  encountered  with  nonconforming  methods  generally  arise 
from  the  nonconformity  of  the  test  space.  Nonconforming  trial  functions  cause 
no  such  problems. 

A  similar  idea  can  be  applied  to  corner  problems.  Consider  solving  Laplace's 

3 

equation  on  a  domain  with  a  comer  angle  of  y  it  ,  and  zero  boundary  conditions. 
The  solution  then  has  a  singular  component  of  the  type  r  sin  2u9/3  .  For 
test  functions  we  use  the  usual  linear  elements,  but  for  trial  functions  we  use 
elements  based  on  the  functions  1  ,  r  sin  2^9/3,  and  r  cos  2^9/3  ,  instead 

of  1  ,  x  ,  and  y  .  Using  this  approach,  the  loss  of  accuracy  due  to  the 
singular  behavior  of  the  solution  in  a  neighborhood  of  the  corner  is  prevented. 

We  remark  that  this  procedure  need  not  entail  any  computational  difficulties 
because  it  can  be  implemented  in  a  way  which  preserves  the  symmetry  of  the  linear 
equations  and  one  may  work  in  the  usual  way  with  the  microstiffness  matrices  and 
nodal  variables. 


.4  CONCLUSIONS 

We  summarize  here  the  main  ideas  we  have  presented.  As  we  have  seen  there 
is  virtually  an  unlimited  variety  of  possible  variational  discretization  methods 
Such  a  method  is  characterized  by  the  bilinear  form  and  the  trial  and  test  space 
In  selecting  a  method  it  is  of  paramount  importance  to  consider  the  goals  of  the 
computation,  in  particular  the  norm  with  which  the  error  is  to  be  assessed.  The 
goals  of  the  computations  are  best  achieved  by  a  method  which  has  good  approxi¬ 
mability  and  stability  properties  with  respect  to  the  desired  norm.  The  method 
should  be  robust  in  the  sense  that  these  properties  apply  uniformly  over  the 


J 


relevant  class  of  problems.  We  note  that  often  the  obvious  method  is  not  best 
and  various  variations  can  lead  to  strikingly  improved  results. 


i i 
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